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Abstract 


A tropical Pacific Ocean General Circulation Model (OGCM) to be used in satellite data 
assimilation studies is described. The transfer of the OGCM from a CYBER-205 at NOAA's 
Geophysical Fluid Dynamics Laboratory to a CRAY-2 at NASA's Ames Research Center is 
documented. Two 3-year model integrations from identical initial conditions but performed on 
those two computers are compared. The model simulations are very similar to each other, as 
expected, but the simulation performed with the higher-precision CRAY-2 is smoother than 
that with the lower-precision CYBER-205. The CYBER-205 and CRAY-2 use 32 and 64-btt 
mantissa arithmetic, respectively. 

The major features of the oceanic circulation in the tropical Pacific, namely the Nort 
Equatorial Current, the North Equatorial Countercurrent, the South Equatonal Current, and the 
Equatorial Undercurrent, are realistically produced and their seasonal cycles are described. The 
OGCM provides a powerful tool for study of tropical oceans and for assimilation of satellite 

altimetry data. 
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1 Introduction 


Over the next decade, many initiatives in remote sensing of the ocean are planned. The 
European Space Agency ERS-1 satellite will be launched in May 1991. The 
TOPEX/POSEIDON altimeter mission is scheduled for launch in the middle of 1992, and more 
altimeter missions are planned in the NASA Earth Observing System (EOS) program. These 
altimeter missions will provide precise and accurate measurements of sea surface height for 

many years. 

Satellite altimeter data, however, only measure information on the ocean surface. 
Techniques have to be developed to extrapolate surface information throughout the ocean 
interior. This is achieved by assimilating satellite measurements into Ocean General Circulation 
Models (OGCMs). The objective of the data assimilation is to obtain dynamically consistent, 
three-dimensional distributions of oceanographic variables throughout the ocean. The resulting 
synoptic description of the ocean is essential for the understanding of ocean circulation and the 
assessment of global climate change. 

Oceanic data assimilation is still in its developing stage and there are few analyses of 
OGCM simulations combined with assimilation of data (Ghil and Malanotte-Rizzoli, 1991). 
Techniques for oceanic data assimilation are mostly confined to simple oceanographic models 
(e.g., shallow-water model or quasi-geostrophic model), which do not adequately describe the 
flow and thermal fields in tropical oceans. Preliminary assimilation experiments of oceanic in 
situ measurements with OGCMs have been started in the Pacific (Derber and Rosati, 1989; 
Leetmaa and Ji, 1989) and in the Atlantic (Morliere et al . , 1989; Carton and Hackert, 1989). 

Satellite data assimilation is currently under way in a joint project at the Jet Propulsion 
Laboratory and University of California, Los Angeles. The project addresses a major question, 
to what extent can the assimilation of satellite sea surface height measurements improve the 

simulation of oceanic circulation by OGCMs? 

The purpose of this report is twofold: (1) to give a brief description of the OGCM, and 

(2) to document the transfer of the OGCM from a CYBER-205 at NOAA’s Geophysical Fluid 
Dynamics Laboratory (GFDL) to a CRAY-2 at NASA’s Ames Research Center. 

2 Description of the Ocean General Circulation Model 

The GFDL OGCM was selected for satellite data assimilation because it can produce 
realistic simulations of current and temperature distributions in the tropical oceans. Over the 
past twenty years, the GFDL OGCM has undergone continuous improvement and is currently 
widely used in the U.S., Europe and Japan. A modified version of the GFDL OGCM 
(Philander et al., 1987) played a very important role in recent developments of dynamic 
oceanography in tropical oceans. This OGCM was first to produce a successful simulation of 
the El Nino Southern Oscillation phenomenon (Philander and Seigel, 1985; Philander et al., 
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1990). A version of this model is operationally run at NOAA's National Meteorological Center 
(NMC); the model output is published each month in the Climate Diagnostic Bulletin, which 
provides a continuing description of monthly oceanographic conditions in the tropical Pacific. 
The GFDL OGCM, in one form or another, is expected to play a central role in the Tropical 
Ocean and Global Atmosphere (TOGA), World Ocean Circulation Experiment (WOCE), and 
EOS programs. 


2.1 Governing equations 


Developments of GFDL OGCM started more than two decades ago. Since then, the 
model has been used to study some of the most basic aspects of large scale, baroclinic oceanic 
circulation. A description of the model physics and numerical methods was published by 
Bryan (1969), and the model code was described by Cox (1984). The model is based on the 
primitive equations with the Boussinessq and hydrostatic approximations. The model's 
prognostic variables are the zonal velocity component (u), meridional velocity component (v), 
temperature (T), and salinity (S). The model's diagnostic variables are vertical velocity 
component (w), pressure (P) and density (p). The governing equations of the model are 
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where A. is longitude, <|> is latitude, z is depth, a is radius of the earth, g is gravitational 
acceleration, Q is angular velocity of the earth, m=sec<]>, and n=sin<|>. The coordinate 
directions K <|> and z are towards east, north and upward, respectively. The three-dimensional 
advection operator is defined as 


T(p) = ma _1 ^(up) + a- I ^(vp) + £-(wp), 


( 2 ) 
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where p. represents u, v, T or S. The horizontal friction terms are 
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where A^ and A TO are the coefficients of horizontal eddy viscosity and diffusivity, and 
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The equation of state is expressed by an empirical polynomial fit to the Knudsen formula 
described in Bryan and Cox (1972). 

2.2 Subgrid-scale parameterizations 

In the governing equations (la) - (lg), the u, v, w, T, and S represent ensemble mean 
variables, which are resolved on the model grid, the primes represent subgrid-scale turbulence 
fields, and overbars represent ensemble means of su bgrid- s cale tu r bulenc e flux es. In order to 
solve the governing equations, turbulent fluxes, i.e., w'u’, w' v', w'T’, and w S , have to 
be parameterized. The first-order turbulence closure scheme assumes that 


w' A' = -K. 


3A 

dz 


(5) 


where A represents u, v, T or S, and K z is the vertical eddy exchange coefficient. In the early 
version of the GFDL OGCM (Bryan 1969; Cox 1984), K z was assumed to be a constant. This 
scheme of constant vertical mixing has problems when it is applied to tropical oceans, because 
the model-simulated ocean is highly sensitive to the K z value (Philander and Pacanowski, 
1980; Semtner and Holland, 1980). For example, although the vertical resolutions in these 
two models are comparable (16 levels in Philander and Pacanowski model and 14 levels in 
Semtner and Holland model), the Equatorial Undercurrent core speed simulated by the constant 
vertical mixing model varies from 0.5 m s 1 in an OGCM with K =1.0x 10' 3 m 2 s 1 (Philander 
and Pacanowski, 1980) to 1.0 m s 1 in an OGCM with K z =1.5x 10 4 m 2 s 1 (Semtner and 
Holland, 1980). The observed core speed of the Equatorial Undercurrent is about 1.0 cm s 
according to Halpem (1987). The main reason for the Equatorial Undercurrent simulated in 
Semtner and Holland model to attain the core speed of 1.0 ms' 1 is the small vertical eddy 
viscosity. 
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Measurements suggested that the vertical eddy viscosity varies considerably in the 
tropical oceans (Halpem, 1980). Given these considerations, a more realistic formulation of 
the first-order turbulence closure scheme should assume that K z depends upon the Richardson 
number, Ri, which takes into account the vertical gradient of temperature and vertical shears of 
currents. The Richardson number is defined as 



where P is the coefficient of the thermal expansion of water and is given as 

P = 8.75 x 10 -6 (T + 9), (7) 

and T is the temperature in degrees Celsius (°C). Pacanowski and Philander (1981) 
formulated the vertical mixing coefficient as an empirical function of Richardson number, given 
as 


K z = K„ + 


Yo 


° (l+oRi)"’ 


( 8 ) 


where K 0 1.34x 10 6 m 2 s y o -5.0x 10' 3 m 2 s' 1 , a=5, and n =2. Given this formulation, 
Philander et al. (1987) showed that tropical oceans can be simulated with a realistic Equatorial 
Undercurrent of 1.0 m s 1 core speed. 

Parameterizations of the constant vertical mixing and Richardson number dependent 
mixing both belong to the first-order turbulence closure scheme. Whatever the empirical 
formula of Richardson number is chosen, the turbulence fluxes are always assumed to be 
locally related to the mean gradient described on the model grid. The subgrid-scale turbulence 
fields, such as the diffusion of turbulence kinetic energy, are neglected. 

Rosati and Miyakoda ( 1988) have adapted a second-order turbulence closure scheme for 
the GFDL OGCM. This scheme, which corresponds to the level 2.5 turbulence closure 
scheme developed by Mellor and Yamada (1973, 1982), includes the turbulence kinetic energy 
as a prognostic variable. Model results with this second-order turbulence closure scheme 
showed a better simulation of the mixed layer in tropical oceans. However, the model 
simulation showed an Equatorial Undercurrent about 50% weaker than that observed. This 
deficiency could be partly due to the vertical resolution in the Rosati and Miyakoda (1988) 
model, where there are only 9 levels in the upper 317 m compared to 18 levels in the Philander 
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et al. (1987) model. Comparison studies of the second-order turbulence closure scheme and 
Richardson number dependent mixing are currently under way (C.-C. Ma, personal 
communication, 1990). 

2.3 Boundary conditions 

At lateral walls, a no-slip condition is applied. No flux of heat and salt is allowed on the 
east and west boundaries. Poleward of 20. 15°S and 32.26°N, the temperature and salinity are 
relaxed to the climatological seasonal values (Levitus, 1982). In these regions, the temperature 
and salinity equations, (lc) and (Id), include an extra term of e(p-p*), where p is the 
climatological value and e is a Newtonian damping coefficient. The values of e depend upon 
latitude and are given in Table 1. This treatment mitigates the effect of artificial zonal walls 
along the southern and northern boundaries. 

At the surface, a rigid lid approximation is made where vertical velocity is zero. The rigid 
lid assumption filters out high speed external gravity waves. The wind stress acts as a body 
force to the first layer of the model. The heat flux (HF) across the air-sea interface is given as 

HF = SW - LW - QS - QL, ( 9 ) 

where SW is short wave radiative heat flux, LW is long wave radiative heat flux, QS is 
sensible heat flux, and QL is latent heat flux. The short wave incoming insolation SW is taken 
to be 242 W m' 1 equatorward of 20° latitude, and to decrease linearly to 145 W m' 1 between 
20° and 45° latitude. The outgoing long wave radiation LW has a constant value of 56 W m l . 
The sensible and latent heat fluxes are calculated using the model- simulated sea surface 
temperature and the prescribed surface air temperature. 

QS = p a C D C P |v|(T 0 -T a ) (10) 


and 


QL = pC D L|v|[E 0 - yE a ](0.622/ 1013.25), ( n > 

where p a is density of air, C D is the drag coefficient and equal to 0.0012, | v| is surface wind 
speed, C p is heat capacity, T 0 is sea surface temperature (°C), T a is surface air temperature 
(°C), L is the latent heat of evaporation and equal to 2.5x 10 6 J kg 1 , and the mixing ratio y is 
fixed as 0.8. The saturation vapor pressure, E, is 

E =10 9 - 4 - 2353/T o, (12) 
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(13) 


E a = 10 9 - 4-2353/1 » . 

2.4 Model design 

The governing equations, along with their boundary equations, are solved numerically by 
finite difference techniques with a staggered "B” grid configuration (Arakawa, 1966), as 
shown in Fig. 1. The time differencing is centered. The model domain covers the Pacific 
Ocean from 30°S to 50°N with realistic coastal geometry, as shown in Fig. 2. The 
longitudinal resolution of the model is 1° longitude. The latitudinal resolution is 1/3° latitude 
within the equatorial band of 10°S and 10°N and increases gradually to 2.5° latitude at 30°S 
and 50°N. The model ocean has constant depth of 4149 m. There are 27 levels in the vertical 
with 18 levels in the upper 317 m (Table 2). The model time step is 1 hour. The CRAY-2 has 
a peak speed of 1.8 billion floating-point operations-per-second (FLOPS). The version of the 
OGCM described herein has an averaged speed of 90 million FLOPS on the 4-processor 
CRAY-2 at Ames. It takes about 6 hours on a CRAY-2 to run one model year. 

3 Comparison of model simulations on CRAY-2 and CYBER-205 

In the numerical experiment to validate the transfer of the model code from GFDL's 
CYBER-205 to Ames s CRAY-2, the Richardson number dependent parameterization is used. 
The initial condition for the experiment consists of climatological-mean January distributions of 
temperature and salinity (Levitus, 1982) with no currents. Climatological-mean monthly 
surface wind stress and air temperature are used to force the model for 3 years. The surface 
wind stress is that of Hellerman and Rosenstein (1983) and the air temperature is that of 
objective analyses of the Comprehensive Ocean-Atmosphere Data Set (Oort et al, 1987). The 
model updates its surface boundary conditions at every time step, which are obtained by linear 
interpolation from the two adjacent monthly means. Snapshots in the middle of each month are 
used in the analyses presented in this section, and 3-day averaged values are used in the next 
section. 

The version of the GFDL OGCM on CYBER-205 used half precision of 64-bit, i.e., 32- 
bit, mantissa arithmetic, which reduced the size of model code by nearly 50% so that it was 
well suited for the 2 million word central memory of the CYBER-205. Because of the 256 
million word central memory of a CRAY-2 supercomputer, the size of the OGCM code is of 
secondary concern so that the full precision of 64-bit mantissa arithmetic is now used. In the 
process of transferring the OGCM code from CYBER-205 to CRAY-2, we performed two 3- 

year model integrations to test the effect of different precision on the simulation of the topical 
Pacific. 

The OGCM code used at Ames is identical to that used at GFDL by Philander et al. 
(1987). The two integrations start from the same initial conditions and are forced with the 
same atmospheric conditions. The differences between the two integrations are compared at 
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the end of month 1, month 13, and month 33 along several selected longitudes (150°E, 180°, 
150°W, 110°W) and latitudes (10°S, EQ, 10°N). Because the same model is used on the 
CRAY-2 and CYBER-205, the differences are due to the accumulation of round-off errors 

introduced by lower precision at every time step. 

The comparison of model simulations on CRAY-2 (represented in the diagrams as A) and 
CYBER-205 (represented as B) at the end of month 1 is shown in Fig. 3 (a), (b) and Fig. 4 
(a), (b) for temperature and zonal current, respectively. Temperature distributions are identical 

and a slight difference is noted in the zonal current distributions. 

The differences between the two simulations are more apparent in month 13, as shown in 
Fig. 5 (a), (b) and Fig. 6 (a), (b). The maximum difference in temperature is located in the 
eastern equatorial Pacific. Differences between the two simulations are larger for zonal current 
than for temperature, with maximum deviation of zonal current occurring at 0° to 5°N along 

150°W. 

In month 33, as shown in Fig. 7 (a) and (b), the maximum difference in temperature is 
between 4°S and 4°N along 150°W. Differences for the zonal current between the two 
simulations are strongest between 3°S to 2°N along 150°W and 110°W (Fig. 8 (a) and (b)). 

The maximum differences of the mean temperature and zonal current along a section 
computed from the two simulations was 0.2°C and 5 cm s" 1 (Tables 3 and 4). Major 
discrepancies occurred in the eastern equatorial Pacific where there is a strong current shear 
between the South Equatorial Current and North Equatorial Countercurrent. The agreement 
between the two simulations was better for month 33 than month 13 along the latitudinal 
sections and also along longitudinal sections of 150°E and 180°, but the agreement was worse 
for month 33 than month 13 along longitudinal sections of 150°W and 1 10°W (Tables 3 and 
4). The average (among seven selected sections) root mean square (rms) difference between 
the temperature signatures was 35% smaller for month 33 than for month 13 (Tables 3 and 4). 
The average rms differences between the zonal current signatures for month 13 and month 33 
were equal to within 5 % (Tables 3 and 4). 

The CRAY-2 simulation has less variability with small spatial scales than the CYBER- 
205 result, although the same model with the same initial and boundary conditions is used. 
The smoother simulated distributions of temperature (Fig. 9) and zonal current (Fig. 10) 
produced by the CRAY-2 compared to the CYBER-205 are caused by the smaller amount of 
computational noise associated with the higher precision of the CRAY-2. Computational noise 
may influence climate predictions (Williamson and Washington, 1973). 

4 Climatological seasonal cycle 

Seasonally, the trade winds along the equator are intense during the northern autumn 
when the Intertropical Convergence Zone (ITCZ) is farthest north, and are weak in northern 
spring when the ITCZ is close to the equator. The westward flowing South Equatorial Current 
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along the equator is strong when the trade wind is intense, and it is weaker when the strength 
of the trade wind is also weaker (Fig. 1 1). During the northern autumn and winter, instability 
waves with oscillation period of 20 to 30 days occur between 0° and 10°N. Signatures of 
these instability waves arc seen in the temperature and current fields. 

The simulated climatological seasonal cycle carried out on the CRAY-2 is similar to the 
CYBER-205 simulation made by Philander et al. (1987). The seasonal cycle of the Equatorial 
Undercurrent is well simulated, as shown in Fig. 12. Similar to observations (Halpem and 
Weisberg, 1989), the Equatorial Undercurrent core is located at about 100 m; its speed is 
strongest (1.2 m s ^) in spring and weakest (0.8 m s"^) in autumn. 

5 Discussion 

A general circulation model of the tropical Pacific Ocean is described. Two 3-year model 
integrations from identical initial conditions but performed on two computers, the CRAY-2 at 
NASA Ames Research Center and the CYBER-205 at GFDL, are compared. The model 
simulations are very similar to each other, as expected, but the simulation performed with the 
higher-precision CRAY-2 is smoother than that with the lower-precision CYBER-205. This 
difference is attributed to the different precision used on the two computing machines. The 
CYBER-205 uses 32-bit mantissa arithmetic and the CRAY-2 uses 64-bit mantissa arithmetic. 
Because of different precision used in the simulations, the round-off errors are accumulated at 
every time step. Comparison of snapshots, on 17th of May and 15th of October during the 
third year simulation, indicates that the differences appeared primarily with small spatial scales 
(Figs. 8 and 9). It appears, therefore, that a significant part of variability with small spatial 
scales simulated by the CYBER-205 is artificially due to truncation errors introduced by the 
lower precision. This computational error may influence climate predictions as detected by 
Williamson and Washington (1973). 

The model described in this report provides a powerful tool for the study of tropical 
oceans, such as the dynamics of the North Equatorial Countercurrent and Equatorial 
Undercurrent. The satellite data assimilation using this model is currently under way. 

At the time of this writing, a new version of the GFDL OGCM (Pacanowski et al, 
1991) has been released. This new model includes essentially the same dynamics described in 
section 2.1. In addition, it includes the three subgrid-scale parameterizations described in 
section 2.2. These different subgrid-scale parameterizations are designed in a modular form so 
that they can be easily applied to various oceanic problems over a wide range of space and time 
scales. The new model code also has a multi-tasking capability which allows about 80% 
parallelism on a four processor CRAY Y-MP. This new version GFDL OGCM will be 
transferred to CRAY-2 at NASA Ames Research Center in the near future. 
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Table 1: Values of Newtonian damping coefficient, e, near the southern and northern 
boundaries. 


j -value 

latitude 

damping coefficient 

(serond) — 

1 

28.84S 

1.98E-6 

2 

26.53S 

2.31E-6 

3 

24.29S 

1.98E-6 

4 

22.15S 

1.16E-6 

5 

20.15S 

3.39E-7 

94 

32.26N 

1.55E-7 

95 

34.78N 

5.79E-7 

96 

37.41N 

1.16E-6 

97 

40.1 IN 

1.74E-6 

98 

42.89N 

2.16E-6 

99 

45.7 IN 

2.31E-6 

100 

48.57N 

2.16E-6 
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Table 2: Vertical levels, grid sizes (m), and corresponding depths of bottom of layer (m). 



12 













































































Table 3: The mean value and root mean square of the difference field between model 
simulations on CRAY-2 and CYBER-205 along 10°S, EQ, and 10°N. 


statistics 

(UCLA-GFDL) 

temperature (degree) 

month 01 

month 13 

month 33 

ION 

mean 

0.00 

-0.01 

0.00 

mis 

0.01 

0.11 

0.01 

EQ 

mean 

0.00 

0.05 

0.06 

rms 

0.02 

0.37 

0.25 

10S 

mean 

0.00 

-0.01 

0.01 

rms 

0.01 

0.11 

0.01 

average 

mean 

0.00 

0.01 

0.02 

rms 

0.01 

0.23 

0.15 


statistics 

zonal current (cm/sec) 

(UCLA-GFDL) 

month 01 

month 13 

month 33 

ION 

mean 

0.67 

-1.19 

-0.03 


rms 

2.15 

4.01 

0.53 

FO 

mean 

0.16 

-3.66 

-0.30 


rms 

3.94 

14.12 

13.19 

ins 

mean 

-0.63 

-1.58 

0.00 


rms 

1.91 

2.22 

0.22 

avnni op 

mean 

0.07 

-2.14 

-0.11 


rms 

2.82 

8.57 

7.62 
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Table 4: The mean value and root mean square of the difference field between model 
simulations on CRAY-2 and CYBER-205 along 150°E, 180°, 150°W, and 110°W. 


statistics 

(UCLA-GFDL) 

temperature (degree) 

month 01 month 13 month 33 

150E 

mean 

0.00 

-0.01 

0.00 

mis 

0.00 

0.09 

0.01 

180 

mean 

0.00 

-0.04 

-0.04 

rms 

0.01 

0.06 

0.05 

150W 

mean 

-0.01 

-0.17 

-0.03 

rms 

0.02 


0.16 

now 

mean 

0.00 

0.13 

0.01 

mis 

0.01 

0.22 

0.11 

average 

mean 

0.00 

-0.02 

-0.02 

rms 

0.01 

0.16 

0.10 


statistics 

(UCLA-GFDL) 

zonal current (cm/sec) 

month 01 

month 13 

month 33 

150E 

mean 

-1.08 

-0.54 

-0.01 

rms 

1.22 

5.86 

0.66 

180 

mean 

-1.24 

-2.01 

0.06 

rms 

4.17 

6.61 

2.14 

150W 

mean 

-3.76 

-3.52 

-4.52 

rms 

4.06 

11.41 

16.82 

now 

mean 

-0.99 

-1.38 

1.87 

rms 

2.18 

5.24 

7.35 

average 

mean 

-1.77 

-1.86 

-0.65 

rms 

3.17 

7.68 

9.25 
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Fig. 2 The model geometry and distributions of model horizontal grids (after Philander et 
ai, 1987). 
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Fig. 3 (a) Surface temperature distributions along 10°S, EQ, and 10°N on 15 February of the 
first year simulated on CRAY-2 at NASA Ames Research Center (represented as A) 
and CYBER- 205 at GFDL (represented as B). 
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Fig. 3 (b) Surface temperature distributions along 150°E, 180°, 150°W, and 110°W on 15 
February of the first year simulated on CRAY-2 at NASA Ames Research Center 
(represented as A) and CYBER-205 at GFDL (represented as B). 






CM/SEC CM/SEC 



LONGITUDE 




Fig. 4 (a) The same as Fig. 3 (a) but for surface zonal current distributions. 
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Fig. 5 (b) The same as Fig. 3 (b) but on 15 February during the second year simulation 
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Fig. 7 (b) The same as Fig. 3 (b) but on 15 October during the third year simulation 
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Fig. 8 (a) The same as Fig. 4 (a) but on 15 October during the third year simulation. 
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Fig. 1 1 Seasonal variations of the surface temperature (°C), zonal and meridional current 
(cm s' 1 ) at 140°W. Positive values in current indicate eastward and northward 
direction. 
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Fig. 12 Seasonal variations of temperature (°C) and zonal current (cm s" 1 ) at 140°W, EQ. 
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